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Several methodologies are developed for large-scale atomistic simulations with fully quantum me- 
chanical description of electron systems. The important methodological concepts are (i) generalized 
Wannier state, (ii) Krylov subspace and (iii) hybrid scheme within quantum mechanics. Test cal- 
culations are done with upto 10^ atoms using a standard workstation. As a practical nanoscale 
calculation, the dynamical fracture of nanocrystalline silicon was simulated. 
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I. INTRODUCTION 

Nanoscale materials are directly governed by quan- 
tum mechanical freedoms of electron systems. Nowa- 
days, electron systems of realistic materials are treated 
by the ah initio electronic structure calculations that are 
based on the density functional theory (DFT)i^, the 
first-principle molecular dynamics'^, and related theories 
developed for decades. A typical system size of present ab 
initio calculations is, however, on the order of 10^ atoms 
and new practical theories are required for nanoscale cal- 
culations. This article is devoted to the methods in large- 
scale electronic structure calculations and their applica- 
tion to nanoscale materials^i^iSi^. 

In general, a quantum mechanical calculation of an 
electron system is reduced to an eigen value equation; 
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with an effective one-body Hamiltonian H. Here the 
eigen energies and eigen states are denoted as {e^'''^^} 

and {0^°'^''}, respectively. A physical quantity (X) is El- We use a standard work station with one Pentium 4^ 
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FIG. 1: The computational time of bulk silicon as the func- 
tion of the number of atoms (N), up to 1,423,909 atoms-; 
The CPU time is measured for one time step in the molecular 
dynamics (MD) simulation. A tight-binding Hamiltonian is 
solved using the exact diagonalization method and an 'order- 
A'^' method with the perturbative Wannier state (See Section 
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with occupied eigen states {(jjf, 
sity matrix p 
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The calculation of eigen states, Eq. is usually reduced 
to a matrix diagonalization procedure and gives a severe 
computational cost. Therefore, the essential methodol- 
ogy for large-scale calculations is how to obtain the den- 
sity matrix p without calculating eigen states. This article 
focuses the methods for structural properties including 
molecular dynamics simulations. For the above purpose, 
the most import physical quantity is the total cnergjiifl. 



The trace in Eq. jSJ is expressed as 

Tr[pX]^ [dr [dr'p{r,r')X{r',r) (4) 
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with real space coordinates r,r'. The off-diagonal com- 
ponents of the density matrix, p{r,r'),r ^ r' , are es- 
sential quantum mechanical freedoms, while the diag- 
onal components, the charge density at the point r 
{p{r,r) = n{r)), appear also in classical mechanics. As 
an important fact for practical large-scale calculations, 
the off-diagonal long range component of the density ma- 
trix dose not contribute explicitly to the value of {X), 
if the operator A" is a short range one. We can found 
a general principle within DFT, called 'nearsightedness 
principle'^^, which is directly related to the above fact. 

For practical algorithms, there are many proposals. 
See reviews or comparison papersi2ii^iiS*iSiii. In this ar- 
ticle, we pick out two methods; (i) method with gener- 
alized Wannier state and (ii) Krylov subspace method. 
Figure ^ demonstrates the computational cost with di- 
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agonalization and our calculation^. Here one can find 
that the diagonalization results in a compuational cost 
proportional to with the system size (N) , as is usual 
in matrix diagonalization procedure (oc N^). Our calcu- 
lation, on the other hand, shows an 'order-A'^' property, 
with upto 10^ atoms, in the sense that the computational 
cost is proportional to the system size (cx N). 

II. GENERALIZED WANNIER STATE 

The generalized Wannier state is a generalization of 
the (conventional) Wannier statesiSti^iSi (See Appendix 
Its pioneering works were done by Walter Kohn in 
the context of large-scale calculationsSiiS^. The pictures 
of the generalized Wannier states are localized 'chemical' 
wave functions in condensed matter, such as a bonding 
orbital or a lone-pair orbital, with a slight spatial exten- 
sion or 'tail'. 

The generalized Wannier states {^l^^"*} are defined as 
localized wave functions that satisfy the equation 

occ 

and the orthogonality 

The matrix Sij is introduced as the Lagrange multiplier 
for the constraint of Eq. © and is given as 

., = (0r)|i/|0r'>. (7) 

The solutions of Eq. (jSJ is equivalent to the unitary trans- 
formation of the eigen states {(l^l^^^H 

occ. 

\c^f^'^)=Y.U^,\4>r'^), (8) 

k 

where Uik is a unitary matrix. Here the suffix i of the 
Wannier state ip^j^^"^ denotes its localization center. 

It is crucial that the generalized Wannier states repro- 
duce the one-body density matrix p in Eq. Q , where the 
eigen states {(/>[f'^^} are replaced by the Wannier states 
{(f)^^^^}- In results, any physical quantity can be repro- 
duced in the trace form of Eq. 101.^^. 

The concept of the generalized Wannier state is used 
for practical large-scale calculationa^i^SiSi. We derived a 
mapped eigen value equation for the generalized Wannier 
states 

^^yr'^)-4^sir^^) (9) 

with a mapped Hamiltonian that is dependent on 

the other Wannier states {(/i^^^^j^^^^. Equation ^ is 
equivalent to Eqs. |SJ| and ©. A Wannier state l^^^^''} 



is not an eigen state of the original Hamiltonian H but 

(i) 

an eigen state of the above mapped Hamiltonian -ff^s- 
Equation ((HJ also shows that the locality of a Wannier 
state can be mapped, formally, to that of a virtual im- 
purity stated. 

With Eq. lO, we developed a variational method so as 
to generate approximate Wannier states, which is called 
variational Wannier state method'*. As the practical pro- 
cedure, Eq. © is solved iteratively under explicit local- 
ization constraint on each Wannier state 

wr'} - {<s} - ir^^} - {<u - ■ • ■ (10) 

With Eq. lO, we also developed a perturbative method 
to generate Wannier states, which is called perturbative 
Wannier state method^. This method corresponds to 
a non- iterative solution of Eq. ©. It is noteworthy 
that the perturbative Wannier state, unlike the varia- 
tional one, is localized without any explicit localization 
constraint, when a short range Hamiltonian H is used. 



III. KRYLOV SUBSPACE METHOD 

The Krylov subspace method^^ is focused as another 
important concept for large-scale calculations. The 
Krylov subspace is a mathematical concept and its def- 
inition is the linear space that is constructed from the 
following vectors; 

\^), H\^), H'\t), ■■■H'^-^\i). (11) 

Here an initial vector \i) should be given. In the present 
context, the matrix iJ is a Hamiltonian. The number of 
bases in the Krylov subspace {v) is chosen to be much 
smaller than that of the original Hamiltonian matrix H . 
In a practical method for large-scale calculations, we 
consider the Hamiltonian matrix only within the above 
subspace, which means the drastic reduction of the ma- 
trix size. The Krylov subspace gives the mathematical 
foundation of many numerical algorithms such as the 
standard conjugate gradient methodSS and the recursion 
method22iffi. 

Recently, we developed a practical Krylov subspace 
method for the calculation of the density matrix and ap- 
plied it to molecular dynamics simulations*. The number 
of bases in the subspace was chosen, typically, as = 30. 
We compared the resultant density matrix with that by 
the Wannier state method. Now we are planning other 
molecular dynamics simulations, especially, to metals. 



IV. HYBRID SCHEME AND PARALLEL 
COMPUTATIONS 

As another fundamental methodology for large-scale 
calculations, we developed the hybrid scheme within 
quantum mechanics^. The basic idea is the foUowings; 




FIG. 2: Snapshots of a fracture process in the (001) plane^. The sample contains 4501 atoms and one initial defect bond as 
the fracture seed. The time interval between two successive snapshots is 0.3 ps, except that between (f) and (g) (approximately 
1.3 ps). A set of connected black rod and black ball corresponds to an asymmetric dimer. See details in the original papej-. 



The one-body density matrix is decomposed into two par- 
tial matrices or 'subsystems' that are constructed from 
several occupied wave functions. This decomposition cor- 
responds to dividing the occupied Hilbert space. The 
different partial density matrices are solved by different 
solver methods. Each subsystem is obtained with a well- 
defined mapped Hamiltonian and a well-defined electron 
number—. Test calculations are done by the combinations 
between (a) the diagonalization method and perturbative 
Wannier state methods, (b) the variational and pertur- 
bative Wannier state methods^ (c) the Krylov subspace 
method and the perturbative Wannier state method^. 
Since the present hybrid scheme is a technique in cal- 
culating the density matrix any physical quantity is 
quantum mechanically well defined with Eq. @. 

Parallel computation is also important for large-scale 
calculations. Test calculation of the perturbative Wan- 
nier state method is carried out with upto 10^ atoms'" us- 
ing the Message Passing Interface techniquei?"'^ and with 
upto 10^ atoms^ using the OpenMP technique^^. We are 
now developing the parallelization of other methodsSiS. 



V. APPLICATION AND DISCUSSION 

As a practical nanoscale application, the molecular dy- 
namics simulation is performed for fracture of nanocrys- 
talline silicon-. A standard workstation is used for the 
simulations with upto 10^ atoms. We use the hybrid 
scheme between the variational Wannier state method 
and the perturbative Wannier state method. 

In the continuum theory of fractura^SiM, a critical 
crack length is defined by a dimensional analysis between 
the competitive energy terms of the bulk strain (3D) en- 
ergy and the surface formation (2D) energy. Since the 
definition of the critical length is independent on the sam- 



ple size or the lattice constant, we can expect a crossover 
in fracture phenomena between nanoscale and macroscale 
samples. The investigation of the above crossover is one 
purpose of the present simulation. Another purpose is 
the fracture behavior in atomistic pictures, on the points 
of how and why the fracture path is formed and propa- 
gates in the crystalline geometry^. 

Dynamical fracture processes are simulated under ex- 
ternal loads in the [001] direction. As the elementary 
process in fracture, we observe a two-stage surface re- 
construction process. The process contains the dras- 
tic change of the Wannier states from the bulk {sp^) 
bonding state to surface ones. Figure |21 shows a re- 
sult, in which the fracture propagates anisotropically on 
the (001) plane and reconstructed surfaces appear with 
asymmetric dimers^. Step structures are formed in larger 
systems so as to reduce the anisotropic surface strain en- 
ergy within a flat (001) surface. Such a step formation 
is understood as the beginning of a crossover between 
nanoscale and macroscale samples^. Further investiga- 
tion should be done for direct discussion of the crossover. 

The present calculations are carried out using tight- 
binding Hamiltonian within s and p orbitals. We should 
say that its applicability is rather limited, due to the sim- 
plicity of Hamiltonian. Its parameter theory, however, 
reproduces systematically several ah initio results among 
different elements or phases, because the tight-binding 
formulation is universal within the scaled length and en- 
ergy units^i^. An important future work is to construct 
simple and practical (tight-binding) Hamiltonians more 
systematically from the ab initio theory. We will use the 
mufhn-tin orbital formulation for the construction, be- 
cause it gives directly the tight-binding formulation^. 

Recently the concept 'multiscale mechanics' is focused 
as the seamless theoretical connection of material simula- 
tion methods among the three principles of mechanics; (I) 
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quantum mechanics (for electron systems), (II) classical 
mechanics, and (III) continuum mechanics. The present 
work gives a guiding principle and a typical example for 
the concept, which is carried out by simplifying the total 
energy functional. 



where the integration is done within the Brillouin zone^. 
Equation (jHJ will be reduced to Eq. l|Aip , when the 
corresponding unitary matrix U is chosen as 



APPENDIX A: DERIVATION OF 
CONVENTIONAL WANNIER STATE 

Here we derive the conventional Wannier stateiSiSS as a 
specific case of Eq. ©. In periodic systems, eigen states 

are called Bloch states {'>pl^°'^^'^ } with the suffices of the 
band v and the k-point k, the point in the Brillouin zone. 
Within an isolated single band, the Wannier states can 
be defined W^i with the suffices of the band ly and the 
lattice vector I, 

W,i{r) = j dke~'^^i^f^°''''\r), (Al) 



U^J U.i.yk = Kv' e-^'^r (A2) 



The conventional Wannier state is given by the unitary 
transform only within an isolated single band {v = i/'), 
while the generalized Wannier states are given by the uni- 
tary transform within different bands {v ^ v'). It is also 
noteworthy that the concept of the generalized Wannier 
state, unlike the conventional one, can be applicable to 
non-periodic cases. 



^ P. Hohenberg and W. Kohn, Phys. Rev. 136, 864 (1964). 
2 W. Kohn and L. S. Sham, Phys. Rev. 140A, 1133 (1965). 
^ R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 

* T. Hoshi and T. Fujiwara, J. Phys. Soc. Jpn. 69, 3773 
(2000). 

^ T. Hoshi and T. Fujiwara, Surf. Sci. 493, 659 (2001). 

® T. Hoshi and T. Fujiwara, preprint (cond-mat 02103661, 

to appear in J. Phys. Soc. Jpn. 72, No.lO (2003) . 

M. Geshi, T. Hoshi, and T. Fujiwara, preprint 

(cond-mat/0306461 1, to appear in J. Phys. Soc. Jpn. 72, 

No.ll (2003) . 

* R. Takayama, T. Hoshi, and T. Fujiwara, in preparation . 
® T. Hoshi, R. Takayama, and T. Fujiwara, in preparation . 

There also exist large-scale calculation methods for trans- 
port properties, whose methodological foundation is quite 
different from the present one. We can find an example^^ 
with the use of the Kubo formula. For review, see the 
following proceedings of a recent international conference; 
RIKEN REVIEW 29, (2000). 

" S. Roche, Phys. Rev. B59, 2284 (1999). 

^2 W. Kohn, Phys. Rev. Lett. 76, 3168 (1996). 

^3 P. Ordejon, Comp. Mat. Sci. 12, 157 (1998). 

" S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999). 
G. GaUi, Phys. Stat. Sol. (b) 217, 231 (2000). 
S. Y. Wu and C. S. Jayanthi, Phys. Rep. 358, 1 (2002). 
D. R. Bowler, M. Aoki, C. M. Goringe, A. P. Horsfield and 
D. G. Pettifor, Modelling Simul. Mater. Sci. Eng. 5, 199 
(1997). 

G. H. Wannier, Phys. Rev. 52, 191 (1937). 
1^ W. Kohn, Phys. Rev. 115, 809 (1959). 

See textbooks, such as N. W. Ashcroft and N. D. Mermin, 
Solid State Physics, Saunders College, Philadelphia (1976). 

21 W. Kohn, Phys. Rev. B7, 4388 (1973). 

22 W. Kohn, Chem. Phys. Lett. 208, 167 (1993). 

2^ In several papers, the Wannier states are constructed in a 
post process of the standard electronic structure calcula- 
tions with eigen states. For silicon crystal, the calculation 
within DFT was carried out using plane wave basea— or 



using muffin-tin orbital bases--^'. Though their methodolo- 
gies in constructing Wannier states are quite different, the 
resultant wave functions show the character of sjf' bonding 
orbital with 'tail'. The 'tail' part of the wave function has 
node structures on neighbor bond sites, as a consequence 
of the orthogonality between the Wannier states. 
2" N. Marzari and D. Vanderbilt, Phys. Rev. B56, 12847 
(1997). 

2^ O. K. Andersen, T. Saha-Dasgupta, and S. Ezhov, BuU. 

Mater. Sci. 26, 19 (2003). 
2** F. Mauri, G. GaUi, and R. Car, Phys. Rev. B47, 9973 

(1993). 

2^ P. Ordejon, D. A. Drabold, M. P. Grumbach, and R. Mar- 
tin, Phys. Rev. B48, 14646 (1993). 

2® As textbooks, G. H. Golub and C. F. Van Loan, Matrix 
Computations, third ed., Johns Hopkins University Press, 
London (1996); H. A. van der Vorst, Iterative Krylov meth- 
ods for large linear systems, Cambridge University Press 
(2003). 

2^ R. Haydock, The recursive solution of the Schrodinger 

equaiton', in Solid state physics, ed. H. Ehrenreich, F. Seitz, 

D. TurnbuU 35, 215 (1980). 
^° D. G. Peffifor, Phys. Rev. Le tt. 63, 2480 (1989). 
^1 jhttp://www.mpi-forum.org/| . 
^2 TittpT/^^'^lopenmp'orgTr"^ 

A. A. Griffith, Philos. Trans. R. Soc. London Ser. A 221, 

163 (1920). 

As a textbook of fracture, see B. Lawn, Fracture of brittle 
solids, 2nd ed., Cambridge University Press (1993) 
See Refi^ and references therein. 

O. K. Andersen, O. Jepsen, and D. Glotzel, in Highlights 

of condensed matter theory. North Holland (1985). 

In general, Bloch states are not unique with respect to 



arbitrary phase freedoms (0^'j!^'' (f") 4>)^^^' (r) exp{i5^k))- 
This phase freedom, however, must be properly chosen to 
construct localized states in Eq. lAll . 



A(Gig) / 



